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Abstract 

A simple but accurate mesh current analysis is performed on a XY 
model and on a SIMF model to derive the equations for a Josephson 
junction array. The equations obtained here turn out to be different 
from other equations already existing in the literature. Moreover, it 
is shown that the two models come from a unique hidden structure. 
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1 Introduction 



This is an analysis of a planar electric network made of Josephson junc- 
tions connected by superconducting wires, biased by external currents in the 
presence of a perpendicular external static magnetic field. Such a circuit is 
usually named Josephson array and it is interesting for technological and the- 
oretical reasons. It is modeled by a system of coupled nonlinear dynamical 
equations, and its physics is based on the behavior of supercurrents interact- 
ing directly each other and also by induced magnetic flux. When there is no 
external bias, starting from any initial condition the supercurrents reach an 
equilibrium state and then circulate forever without dissipation, sometimes 
trapping flux in the network meshes and accomodating themselves in char- 
acteristic patterns which also depend on the external magnetic field. In the 
presence of bias these currents can rearrange and adjust to the new external 
conditions, up to a certain bias limit. If this limit is exceeded by the bias, 
an out of equilibrium dissipative state is reached and flux patterns can move 
through. The circuit then emits a radiation originated in its elements: if syn- 
chronization sets up, this radiation can be constructively coherent and much 
more energetic than the radiation from a single junction. The exploitation 
of this synchronization mechanism is the technical reason for the interest for 
such a system, which allows making usable the feeble but very high frequency 
and easily tunable radiation emitted from a single Josephson junction. 

The theoretical interest is in the understanding of the complex dynamics 
underlying the dance of the flux patterns and the action of the magnetic field 
on both static and dynamic state. To achieve this result, the attention was 
focused on the choice of the model to be used to describe this system. As 
pointed out in Ref.|l|, two seemingly different classes of models have been 
selected: the first one is referred to as XY model and the second one can be 
referred to as self-induced magnetic field (SIMF) model. The XY model de- 
rives from field theory and indirectly explains many array phenomena, but 
is limited to systems in wich array inductance is negligible. The SIMF model 
takes into account the induced flux interaction and then should be valid in 
every situation, including the case when the inductance goes smoothly to 
zero. Exactly at zero inductance it should merge with the XY model, but 
this cannot be proved because of the presence of a parameter singularity. 
The SIMF model is more difficult to be studied analytically than the XY 
model, but in the form studied in Ref.[§] it has a rare direct experimental 
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evidence [3,0. It is anyway thought to be fundamental to explain some 
peculiar features in giant Shapiro steps and other related phenomena. In 
Ref. 1^ a compact formalization of the SIMF model was given as an explicit 
dynamical system. In Ref.|^ an attempt was made to study analitically the 
parameter region where the two models should merge, but the singularity 
in the inductance parameter makes the analysis difficult. A numerical in- 
vestigation of this parameter region was made in Ref. 0, but a numerical 
simulation drawback is the difficulty to get a global understanding of how 
the transition works, beeing limited only to finite parameter ranges. Ref.p| 
is also a source for more references. 

In this work two main conclusions about these models are reached. A 
simple but careful derivation of the equations leads to systems which are 
different from those existing in the previous literature. Moreover, XY and 
SIMF models derive from the same set of equations and the study of the 
transition region is easier than it was thought before. The plan of the paper 
is to work out a simple reference example in Section 2 and to derive in a 
general way the equations for a XY and a SIMF rectangular network in 
Section 3. Then in Section 4 the source of the difference will be briefly 
discussed. In Section 5 the system from which both models derive will be 
shown. 

2 Four circuits as an example 

In principle, a Josephson array can be built laying a rectangular grid of 
superconducting stripes on an insulating substrate so that cells surrounded 
by superconducting wires are formed. A rectangular array has A^,. x Nc cells, 
where Nr indicates the number of cells in one row and Nc the number of 
cells in one column. In the examples of Fig.|l| the cells are square and each 
circuit consists in two cells. In circuit a) conventions are displayed: nodes 
are labeled by the two row and column indices (r, c) and reference directions 
for the currents are chosen. Orientated currents will be called iv,r,c and 
orientated gauge invariant phases will be called (pv^r,c , where r and c refer to 
the node where they originate from, index f = is for horizontal branches 
and V = 1 for vertical branches. Meshes are labeled only with the two indices 
(r, c) of their bottom left hand corner and orientated clockwise. External bias 
currents are called Ir,c where the indices r and c refer to the node where the 
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Figure 1: Example circuits and notation 



current is injected. The Josephson junctions considered here are proximity 
effect junctions, which can be obtained cutting the branch wire in the center, 
and they are indicated by crosses. Small boxes will be explained later. The 
current i flowing along a branch containing a junction is described by the 
resistively shunted (RSJ) model as 
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in terms of the superconductor order parameter phase difference ip between 
the edges of the junction: i? is a phenomenological resistence and ic is the 
maximum supercurrent that can flow in the junction. The first z„ term arises 



from the voltage V 
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if which sets in due to a.c. Josephson effect when 



if varies in time, the second ig term describes the V = Q supercurrent that 
flows by Cooper pairs tunneling. To treat as a real variable and not as a 
distribution, an assumption of very low temperature conditions is made here. 
Three sets of equations will now be written |^ for each of the circuits 6), c), 
d) and e) in Fig.[l|. Since circuits 6), d) and e) have junctions in every branch 
they will be named pure Josephson networks while circuit c) will be named 
hybrid, lacking junctions in some branches. 

When currents circulate in an array, each wire induces a magnetic field 
proportional to the current i flowing in the wire itself. If this magnetic field 
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is not too strong, each cell is threaded only by the self-induced magnetic flux 
(j)r,c generated by the closed current formed by the currents i^,r,c circulating 
in the wires along the borders of the mesh (r, c) , plus a possible external 
magnetic flux 0^^*. If 0*°* is the total flux, 



I tot //)^^* 



r,c (2) 

For mesh (1, 1) and (1,2) in circuit b) induced fluxes are 

f 01,1 = -L(ii,i,i + Zo,2,l - «1,1,2 - «0,l,l) /g\ 
I 01,2 = iv(ii,i,2 + io,2,2 - H,l,3 - ^0,1,2) 

where L is the cell self- inductance. In the circuits in Fig.^some branches have 
small boxes attached to indicate generation of flux, while branches without 
boxes generate negligible flux. If a junction is present in a wire, it limits to 
ic the current that can flow in it; if the junction critical current ic is very 
low, the wire contributes to induced flux only with the z„ normal term, which 
at equilibrium is null. In fact, junctions not enclosed in boxes have critical 
currents i'^ negligible in comparison to others. 

In this aspect hybrid array c) and pure array d) are analogous: in c) 
vertical branches contributions to the flux are negligible in comparison to 
contributions from unbroken horizontal wires. In d) junctions with very low 
i'^ < ic are present in vertical branches and their induced flux is negligible in 
comparison to that from horizontal branches. Then for circuits c) and d) 

f 01,1 = -^(^0,2,1 - «0,l,l) /^N 
\ 01,2 = L{io^2,2 - ^0,1,2) 

while in circuit e) 

01,1 = 01,2 = (5) 

which could be naively seen as a L ^ limit of Eq.(|]). Eq.(|D, (^) and (|^) 
are different instances of the first set of equations. Circuits b), c) and d) are 
SIMP models while circuit e) is a prototypical XY model. 

The second set of equations comes from the principle that the directed 
sum of gauge invariant phase differences around a mesh is related to the 
induced flux threading the mesh 0: 

j 01,1 = -(V^l,l,l + ^^0,2,1 - ^^1,1,2 - V^O,l,l) + /l,l /pN 
1 01,2 = -{^1,1,2 + ^0,2,2 - V'1,1,3 - ^0,1,2) + fl,2 
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where fr^c = 27r rir^c + 0r^* are called for historical reasons frustrations and 
nr,c are integers. In ,0i,2, 4>Tc ^ coefficient ^ has been absorbed, where 
e is the electron charge and h the Plank constant. The sums {<f 1,1,1 + f 0,2,1 — 

^1,1,2 - '^0,1,1 + 01,1 - ^r^c ) (^1,1,2 + ^0,2,2 - ^1,1,3 - ^0,1,2 + 01,2 -0rfc*) 

are named fluxoids and Eq.(|]) states that fluxoids are equal to 2tt times an 
arbitrary integer n. 

The third set of equations implements Kirchhoff current conservation law 
for nodes in an electric network. The oriented currents sum to zero for each 
node as 



which is valid for every circuit. 

From these three sets of equations a nonlinear differential system for 
each circuit is going to be derived, where only one phase derivative will 
appear explicitly for each equation. Such a system is called explicit dynamical 
system. For simplicity, from now until differently stated, it will be assumed 
that external fluxes 0^^* and integers nr,c are zero, and R = 1. A bias current 
7 enters from the upper nodes and leaves the circuit from the lower nodes 
as -7: then l2,i = 12,2 = ^2,3 = 7 and Ii,i = Ii,2 = /i,3 = -7- System (0) 
consists of 6 linear equations in 7 unknowns, but only 5 equations are linearly 
independent. Two free parameters and I12 are introduced, named mesh 
currents, to solve it: 



Direct substitution in system works as a proof. Induced fluxes in systems 
(H) and (|) can be expressed in terms of Eq.(||) too: 



An immediate remark has to be made here: the system from @ relates the 
self-induced flux (j)r,c with all the possible mesh currents whereas the system 




(7) 






(9) 
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from relates the flux (j)r,c only to the mesh current I^^ bearing the same 
indices (r, c) . These systems can be inverted to give 





1 1 



trj'' (10) 



and such mesh currents can be inserted in Eq.(D, which now relates branch 
currents iv^r,c to self-induced fluxes (j)r,c- 

The differential equations for the circuit b) can now be written substitut- 
ing in system (H) Eq.(0) for the l.h.s. and the first block of Eq.(|l^) for the 
r.h.s., where 0r,c are obtained in turn from Eq.(|^ 

V'o,!,! +icSin(v9o,i,i) = 

1^(4V20,1,1 + V50,l,2 - 4V?1,1,1 + 3^1,1,2 + <^l,l,3 - 4V?0,2,1 - <^0,2,2) 
<^0,1,2 +ic Sin((y90,l,2) = 

1^(V^0,1,1 + 4990,1,2 - - V^l,l,2 + 4(y9i,i,3 - (y9o,2,l - 4^90,2,2) 

'/'l,!,! +icSin((y9i,i,i) = 

igl (-4^50,1,1 - V50,l,2 + 4V51,1,1 - 3(/3i_i,2 - V5l,l,3 + 4V90,2,1 + V50,2,2) " 7 

¥^1,1,2 +icsin(v?i,i,2) = 

i^(39?o,i,i - 3v9o,i,2 - 3v2i,i,i + 6v?i,i,2 - 3v?i,i,3 - 3v?o,2,i + 3920,2,2) - 7 
<y^i,i,3 +icsin(99i,i,3) = 

1^(¥'0,1,1 + 4V90,1,2 - V^l,!,! - 3V31,1,2 + 4v2i,i,3 - V5o,2,l - 4V50,2,2) " 7 

V'0,2,1 +icSin(v9o,2,i) = 

pi:(-4vi'o,i,i - V^o,i,2 + 4v2i,i,i - 3v9i,i,2 - </3i,i,3 + 4v?o,2,i + ^^0,2,2) 

V'0,2,2 +«cSin(v9o,2,2) = 

^^(-V^O,!,! - 4V50,1,2 + V^l,l,l + 3V51,1,2 - 4v9i,i,3 + V9o,2,l + 4v9o,2,2) 

This result is remarkably different in the form of the interaction among phases 
from those in Ref.0 and [0. Here in each equation are present a// the phases 
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of the system, while there in each equation is present only a subset of the 
whole set of phases. 

The differential equations for the hybrid circuit c) are obtained in a 
slightly different way, because orizontal branch currents io,r,c are not de- 
scribed by Eq.(0). In this case branches are considered as containing induc- 
tances and not as simple shorts. Substituting 0j.,c from Eq.(|D in the second 
block of Eq.(0) gives /"^ in terms of i^,r,c- In turn are inserted in system 
(g). A set of 4 relations 2o,i,i = -|(^o,2,i - ^0,2,1 = |(^o,2,i - «o,i,i), 
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c) 



io,i,2 = -|(^o,2,2 - "^0,1,2) and ^0,2,2 = |(«o,2,2 - ^0,1,2) is obtained and it is sat- 
isfied iff zo,i,i = ""^0,2,1 and io,i,2 = — "^0,2,2- This implies that the 6 equations 
in system (J^) are reduced to the 3 equations 

7 

(11) 

where Eq.(0) is to be apphed only to the l.h.s. terms. Eq.(||) for circuit c) 
has to be modified in = -{f 1,1,1 - '^1,1,2) and 0i,2 = -{f 1,1,2 - <^i,i,3), 
because only two junctions are present in each cell. Using this relation gives 

V'l,!,! +icsin(v9i,i,i) = -^(V5i,i,i - f 1,1,2) - 7 
'^1,1,2 +icsin(v9m,2) = ^{^1,1,1 - 2v5i,i,2 + V^i.i.s) - 7 
. ^1,1,3 +icsin(v9i,i,3) = ji{^i,i,2 - ^1,1,3) - 7 

This result is the same as that in Ref . , where it was derived for an infinite 
stripe of cells. For this circuit, in contrast to the circuit 6), in each equation 
is present only a subset of all possible phases. This is clearly due to the fact 
that in every mesh of this circuit the self-induced flux 0r,c is proportional 
only to the I^^ with the same indices (r, c). 

The differential equations for the pure array d) are derived in the same 
way as for circuit 6), substituting in system (H) Eq.(|l]) for the l.h.s. and the 
second block of Eq.(|TU|) for the r.h.s., where in turn are obtained from 
Eq.(l): 

</'o,i,i +«cSin(v9o,i,i) = 

~J,\ (V'O,!,! - V'l,!,! + V'1,1,2 - </2o,2,l) 

Vo,i,2 +ic sin(v?o,i,2) = 

-t\ {^0,1,2 - ^1,1,2 + 1,1,3 - ^0,2,2) 



d) 



+icSin(v9i,i,i) 



j_ 1 

' L 2 



(-V^o,i,i + ~ V^i,i,2 + V^o,2,i) - 7 



Vl,l,2 +i'^siYl{ipi,i,2) = 

~Z\ (V^0,1,1 ~ V^0,l,2 - V'l,!,! + 2V51,1,2 - V^l,l,3 " ¥^0,2,1 + V'0,2,2) " 7 

V'l.i.s +^cSin(v?i,i,3) = 

~t\ (v'0,1,2 - V'1,1,2 + V^i,i,3 - V'0,2,2) - 7 

V'0,2,1 +icSin(v?o,2,i) = 

~T.\ (^V^o,i,i + V^i,i,i ~ Vl,l,2 + </50,2,l) 

V'0,2,2 +ic sin(<yCo,2,2) = 



■J^\ [-^0,1,2 + V'1,1,2 - Vl,l,3 + V'0,2,2) 
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This system has the same polarization scheme as systems b) and c) and a 
hmited interaction among phases: only in the differential equation for phase 
V?i,i,2 all the other phases appear. This is again due to the proportionality of 
mesh currents and induced fluxes. It seems having never been discussed in 
the literature, even though it has an interaction equal to that attributed to 
circuit b) in Ref. 0. 

The differential equations for the last circuit e) are obtained in another 
different way. In this case induced fluxes are zero, so it is not possible to use 
them to eliminate the mesh currents in system (|^). Anyway, Eq.(^) is still 
valid and can be differentiated. Inserting Eq.(|l]) in system (|[) together with 
differentiated Eq.(|D gives a system that can be rewritten as 



V'o,!,! +icSin(v9o,i,i) 
'f 0,1,2 +icSin(v?o,i,2) 
V'l,!,! +icsin(v9i,i,i) 
'f 1,1,2 +icsin(v9i,i,2) 
V'l.i.s +icsin(v9i,i,3) 



1,1 

_ 

1,2 

-'1,1 7 

-'1,1 ^ -'1,2 
-'1,2 7 
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'f 0,2,1 + ic sm{(fo,2,i) = III 
V^o,2,2 + ic sin(v9o,2,2) = /l,2 

^^1,1,1 + '^0,2,1 - 7^1,1,2 

- y^o,i,i= 

¥^1,1,2 + 0,2,2 - V^l,l,3 

- ¥^0,1,2= 



Eliminating the last two equations allows to obtain explicitly the mesh cur- 
rents that are now expressed in terms of sine functions of phases: 

Ii,i = ft(-4sin(v9o,i,i) - sin(v9o,i,2) + 4 sin(v9i,i,i) + 
-3sin(v9i,i,2) - sin(v9i,i,3) + 4 sin(v9o,2,i) - sin(v9o,2,2)) 

-^1,2 = ft(-sin(v9o,i,i) - 4sin(v?o,i,2) + sin(v9i_i,i) + 
+ sin(v9i,i,2) - 4 sin(v9i,i,3) + sin(v9o,2,i) + 4 sin(v?o,2,2)) 

After inserting them in Eq(D, an explicit system of 7 differential equations 
(of which only 5 are independent) results. To give a synthetic idea of its 
structure, only two equations are presented here : 

^0,1,1 +icSin(v9o,i,i) = Y|(4sin(v9o,i,i) + sin(v9o,i,2) - 4 sin(v9i,i,i) + 
+3sin(v9i,i,2) + sin(v9i,i,3) - 4 sin(v9o,2,i) - sin(v9o,2,2)) 

^^0,2,2 +icSin(v9o,2,2) = f| (- sin(v9o,i,i) - 4 sin(v?o,i,2) + sin(v9i,i,i) + 
+3 sin(v9i,i,2) - 4 sin(v9i,i,3) + sin(v?o,2,i) + 4 sin(v9o,2,2)) 

Remarkably, this system is very different from the previous ones, having no 
term linear in v^u.r.c- It has the same polarization scheme as the other circuits 
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and the same coefficients as b) in front of the interaction terms, which are 
now sine functions of phases instead of beeing simply phases. It is definitely 
not obtained as a direct L — > limit from system b) . 

3 General derivation of the equations 

One of the two main points of this paper is that a careful derivation of the 
equations for the circuits b), d) and e) leads to systems of equations different 
from those already existing in the literature. In the previous example it 
was shown that a series of manipulations leads to an interaction that has a 
range longer than what precedently thought. This point will be discussed 
more explicitly in the next section. In Section 5 it will be shown that, even 
though the three circuits seem to be rather different systems, their equations 
can be derived from a unique hidden structure and that it is misleading to 
consider the limit L — > directly from equations like those just obtained. To 
discuss more easily the physics involved, two steps will be followed. First, 
the generalization of the preceding derivation will be illustrated, then a new 
set of variables will be discussed. In these variables it will be possible to 
smoothly handle the limit L — > 0. 
For any x Nc array the relation 

n^,- Un+l (12) 

holds, where = N^Nc is the number of meshes, rib = 2NrNc + Nr + 
is the number of branches and n„ = (A^,. + l){Nf. + 1) is the number of 
nodes. Using a matrix notation, the current vector i = {iv,r,c} = {ik, k = 
l,...,nb}, the phase vector = {'fv,r,c} = {fk,k = l,...,nb}, the mesh 
current vector 7" = {-^r,c} = {^k^ k = 1, . . . ,71^}, the external bias current 
vector Ta = {Ir,c} — i^k^k — l,...,nm} and the frustration vector / = 
{/r,c} = {/fcj k = 1, . . . jfim} are introduced . In current and phase vectors 
the elements are progressively ordered by an index k obtained from the indices 
of their elements as A; = {2Nc + l)(r — 1) + vNc + c, r = 1, . . . , iV^ + (1 — v), 
c — 1, . . . , Nr + v: for example, current ii,i,2 of circuit a) becomes ^4. In 
words, starting from the bottom left hand corner and following the row of 
nodes, first come the horizontal branches then come the vertical branches, 
and then next node row is scanned. In frustration and external current 
vectors the elements are progressively ordered by an index k = Nc{r — 1) +c. 
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where r = 1 . . . , Nr and c = 1, . . . , Nc, and the same order is assigned to the 
meshes, considered as rings clockwise oriented. Two matrices Aa and Ma 
are now built |TD|. The matrix Aa, with n„ rows and columns, describes 
the Kirchhoff law for current conservation at the nodes. Its element {Aa)m,n 
is 1 if in the m-th node the n-th current enters the node, it is —1 if the 
current leaves it, it is otherwise. The matrix Ma, with + 1 rows and 
rih columns, describes the sums around meshes. Its element {Ma)m,n is 1 
if in the m-th mesh the n-th current belongs to the mesh and is oriented 
in the same direction of the mesh, it is —1 if the current is oriented in the 
opposite direction and otherwise. The matrix Ma has one row more than 
Um because it includes the so called external mesh, which follows the rule 
that {Ma)ni,+i,n IS —1 if the n-th current is on the boundary of the array and 
is oriented clockwise following the boundary ring, 1 if it is oriented in the 
opposite direction and otherwise. It can be proved using network theory 
[|T^ that rows of Aa and Ma are mutually ortogonal, and this is also seen by 



direct inspection. Since rank(y4a) = n„ — 1 and rank(Ma) = rim, from Eq. (|T^) 
rank(74a)+rank(Ma) = n^ follows, which means that Aa and Ma divide the 
Ub dimensional space in two parts. Deleting a row in each of the two matrices 
does not change their rank: as a convention the last row will be deleted in 
Am to form a new matrix A, and the last row will be deleted in Mm to form 
a new matrix M. These two matrices have maximum rank and the property 

/ MA^ = 0, AM^ = , . 

\ M^{MM^)-^M + A^{AA^)~-^A = 1 ^^'^^ 

having AA^ and MM'^ non-zero determinant. They appear to be dual in 
the sense of linear algebra and the quantities 



M^{MM^)-^M = K 
A^{AA^)-^A = K 



TiAAT^-iA_i7 (14) 



can be seen as projectors on the n^ dimensional space. The induced flux 
vector can be introduced as 

= M^i (15) 

where the self-induced flux matrix M^ expresses a relation analogous to 
Eq.(|^) or Eq.(|^). This matrix is built from M substituting in M for ev- 
ery non-zero entry the branch inductance Xm,n if the entry is 1 and — A,^ 
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if the entry is —1. Eq.(^) is recovered when all Xm,n are equal to L, that 
is = LM. In the general case M^A^ ^ and A{M^Y ^ 0. The last 
element is deleted from vector to form vector T. Since elements obey 
current conservation law -^fc = 0, the deleted element is always obtained 
as = —YJkZi^^k- In this notation equations Eq.(0), Eq.(|^) and Eq.(l^) 
can be rewritten as 

Ai = T 

M^i + Mip = f (16) 

where in sinfc(</))fc the label k is to remind the vectorial form. It is stressed 
that notation Fsmk{G)k indicates that the generic row k of matrix G ap- 
pears as the argument of the k-th element of a vector defined as sin k{G)k = 

(J^ '-j^ 1_ ]_ 2 
J^2,l -^2,2 

An ( i^i. TP ■ (n\ I ^1,1 si'^ '^i + ^1-2 sin \ ^, , 

and G = „ then Fsmfc(G)fc = . Clearly, 



^ -^2,1 sin Gi + F2,2 sin G2 
k is not an index. The equations for circuit h) are now quickly rederived. 
The general solution oi Ai = T is the sum of the solution of the homoge- 
neous equation Ax = and a particular solution of Ax = T. Using Eq.(|T3|), 
i = M'^r + A^D where D = {AA'^y^T. Operating with on both sides 
of this equation resuhs in M^i = M^M'^P + M^A^D, Eq.®, that after in- 
version gives = {M^M^)-^M^i-{M^M^)-^M^A^D, Eq.(^. From the 
second equation of system ( |16|) it follows that M^i = f — Mip, Eq.(||), which 
inserted in the equation for I" gives I" = -{M^M'^)-^Mip + [M^M'^y^f - 
{M^M'^)~^M^A^D. Finally, using the Josephson equation and defining 

K = M^{MM^)-^M 

= M^{M^M^)-^M (17) 

the equation 

+ic sin fc((^)fe = -K^if + M^{M^M^)-^f - {K^^ - l)A^D (18) 

for systems h) and d) is recovered. A comment about the conditions under 
which M^M^ can be inverted is needed. As a general rule det (M^M"^) 7^ 
if at least rim branch inductances are non-zero and every mesh has at least 
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one branch with non-zero inductance. This is the case of circuits b) and d) 
but not of circuit e). Moreover, it can be shown that for any array dimension 
no element of K is zero, so that K generates an infinite range interaction 
among phases. 

In the case of circuit e) starting equations are 

' Ai = T 

< Mif = f (19) 

If / = it follows that M ^= 0. Operating with M on the Joseph- 
son equation gives M M icsin k{if)k + Mi, where the l.h.s. member 
is zero. Inserting in this equation the solution i = M'^I"' + A^D gives 
MicSin k{(p)k = MM'^I"', where Eq.(0) has been used. After inversion, this 
results in I"^ = {MM'^)~^MicSm k{(p)k- Using again the Josephson equa- 
tion +ic sin = M'^I" + A^D and Eq. ([T3|) the equation for e) quickly 
follows: 

^ +icsink{(f)k = icK sin k{ip)k + A^ D (20) 

If / 7^ this method doesn't work, because it gives again the same equation 
without any trace of /. The use of Myj = / as a constraint is critical. When 
differentiated this equation is always M V?= independently of /. This 
problem can be avoided introducing another variable p as already done with 
This variable will be called cut phase, in analogy with network theory 
terminology. More precisely, the equation Myj = / can be solved by the 
substitution = A^p + M'^ s where s = {MM^)~^ f . The quantity M is 
still zero but at the end of the calculation the result is 

A'^ P= {K -l)i^sink{A^p + M^{MM^)-^f)k + A^D (21) 

which can be multiplied by A and inverted to give 

P= -{AA^Y^A i,sink{A^p + M^{MM^)-^f)k + {AA^y^T (22) 

where Eq. ([13|) has been used. After solving this differential system, phases ip 
are recovered as ip{t) = A^piJ:) + M'^s. This new variable p seems to better 
suit to the L = problem, transforming it in a system of independent and 
explicit differential equations. Unlike the case of Eq. ([T6D it is not possible to 
eliminate p at the end of the calculation, as it was did with the intermediate 
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variable The reason stands on the fact that system (|T9|) is constrained 
by the rim static equations Mip = f and only nt — rim variables are left inde- 
pendent, so that the whole dynamics depends only on these new variables. 
Applying this new point of view to the circuit e) , the system of 5 independent 
explicit differential equations in the 5 (/ = 1, . . . , 5) variables 



(Pi \ 

Ps 

Pi 
\ P5 ) 



^ sin( 
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-5 
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10 


10 


1 


4 


-1 


-3 


4 


-1 


11 



/ 7 \ 

7 
7 


(23) 
results. 

This idea suggests a way to write the general system in a form suitable to 
take the L — > limit in a smooth way, but before undertaking this task the 
difference between Eq.(|T8l) and other existing equations will be discussed. 



sm(^-pi +P2) \ 

sin(-p2 +P3) 
sin(-pi +P4) 

sin(-p2 +P5) 

sin(-p3) 
sin(-p4 +P5) 
sin(-p5) ) 



4 Mesh currents and induced fluxes 



It is stated in Eq.(6) from Ref.|[TT| that for an array where only the self- 
inductance of each lattice mesh is retained, and the mutual inductance among 
cells is ignored, the relation between the total flux in one mesh and the 
mesh current of the same mesh Jjj is 



itot 



t] r'ext " ^ij 



L'l-, (24) 



This means that the self-induced flux is proportional to the mesh current. 
This equation is used also in Ref. [0 and |T^ and implicitly in Ref . |T3[ . In 



the present work it has been shown for hybrid circuit c), which is a particular 
1x2 case of the 1 x 00 array of 0, that Eq.(^) is correct because in 
this particular case the induced fluxes are indeed proportional to the mesh 
currents. From the second block of Eq. ([To|) it is seen that Eq.(p4|) is also valid 
for the pure circuit d), but it is not in general valid, e.g. when the magnetic 
fields induced by the four branch currents surrounding a cell are all of the 
same strenght. This happens for example in circuit b), as it is seen from 
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the first block of Eg. (p!OD . This is the source of the difference between the 
equation for circuit b) derived here and the equations derived using Eq.(2^). 
Yet, tliere is a case in wliicli Eq.(^) can be used also for circuits similar to 
b) where all the four branches contribute in the same way to the flux. Now 
it will be shown that in the formalism developed in this paper it is easy to 
work out an interaction in which Eq.(^) is anyway valid. To achieve this 
goal, the mutual inductance among meshes has to be taken into account. 

The task is to link mesh currents and mutually induced fluxes in such a 
way that 

(Pr,c = M^'i = M^'M^r = p I^^^ (25) 

where in analogy with Eq.(|T5|) M^^^ is a mutually induced flux matrix and p 
is an adjustable constant. For circuit b) in the example of Fig.|l|, 



M 



MI 



A^^ 



XFN 



where A"^-^ is a self-inductance and A^^ is a first neighbour mutual inductance, 
and 

' -1 1-1 1 \ 
-10 1 -1 1 J 



M 



Requiring that M^'^ 
system 



pi, where 1 is the unitary matrix, gives the 



4^57 _ ^FN 
^^FN _ ^SI 

4^57 _ 



1 



1 



A 



SI 



When such an interaction is 



15 

4 ^ 4 ■ 

is recovered. This model can be called MIMF 



and p 



which is solved by A'^^ 
taken into account, Eq.l 
(Mutually Induced Magnetic Flux) model. From this result an immediate 
consideration follows. In a circuit like b) with only self- inductance taken into 
account, the interaction among phases is global: every phase of the system 
appears in each differential equation. It is interesting to realize that a local 
phase interaction typical of a circuit like d), where in each equation only 
a small subset of the phases appears, can be obtained introducing in b) a 
global mutual inductance coupling. In this sense interaction range in mutual 
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inductance among cells and interaction range in coupling among phases are 
dual. 

5 The L ^ limit 

To focus attention on the structure of the limit, it will be assumed that 
in Eq.(]T6|) = LM, which means that circuit b) will now be studied. 
Equation M[Li + (p) = f can be solved by setting q = Li + (f where q = 
A^p + M'^s. Equation Ai = T is solved hj i = 1°" + Al^ and merging 
these two relations gives 

^ = A^p - LM^I" + M^s - LA^D (26) 

Inserting this definition in the Josephson relation gives 

A^ P -LM^ I" +ic sin kiA^p - LM^P + M^s - LA^D)k = M^P + A^D 

(27) 

In a block matrix form, defining the matrices 

= (A^, -LM^), B'^ = (0, M^) (28) 

and the vector c = {p, P), Eg . (|27|) is rewritten as 

B^ c= -ic sin fc(5^c + M^s - LA^D)k + B^\ + A^ D (29) 

which after inversion is the differential system 

P= - {AA^)-^A sin fc(A^p - LM^P + M^s - LA^D)k + D 
L f= (MM^)-^Mie sin fc(A^p - LM^ I" + M^s - LA^D)k - I" 
D = {AA^y^T, s = {MM^)-^f 

(30) 

In these variables the limit L — is easier to handle than in a form like 

f+i^sm{^) = ^^ (31) 

In the first row of Eq.(^) the interaction between p and LP smoothly van- 
ishes and at L = Eq.(p2|) is recovered. In the second row the dynamics 
gets smoothly frozen and at L = a static constraint P = {MM^)~^icM 
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sin k{A'^ p + M'^s)k is imposed to the space of the degrees of freedom (d.o.f.). 
The variables p dynamically decouple from and depend statically on 
p. Multiplying this equations by MM'^ gives Mi = sin ^(A-^p + M^s)^ 
which going back to ip is the condition M{i — ic sin{(f)) = that follows from 
M 'P= 0. The rib dimensional d.o.f. space is separated in two sectors, the 
dimensional mesh currents sector and the Uf, — rim dimensional cut phases 
sector, and dynamics survives only in the cut phases sector. This is a nice 
example of what could be called a parametric dynamical bifurcation, in anal- 
ogy with the usual parametric Hopf-like bifurcation of dynamical system 



theory. When this differential system is to be studied by perturbation theory 
a singular problem is obtained, but in this form the stiffness is less harm- 
ful. When L is very small but not zero, the time variable t in the derivative 
d/dt can be replaced by t/L to use a dilated time scale. In this scale the 
p appear almost frozen. The almost slaved tend to follow the dynamics of 
the p, yet display a second faster but strongly damped dynamics. Moreover, 
the behaviour of Eq. (|30|) under L — > limit is physically satisfactory, because 
it is hardly belivable that changing smoothly the inductance of these circuits, 
e.g. making them bigger, should give some sharp effect only and exactly at 
L = 0. When instead of LM is used, notation becomes heavier, but it 
can be shown that only belonging to meshes with zero inductances around 
the cell become constraints while the remaining keep on taking part in 
the dynamics. 

Eq. (^) hides a nice asymmetry or duality between external currents and 
external magnetic fluxes. To make it explicit, for L 7^ the change of 
variables is made 

(u=p-LD , . 

\w = -LI- + s 

Multiplying the first row of Eq. ( |30D by and the second row by M^, and 
substituting in it Eq.(^), gives 

J u= -Kic sin ^(A^m + LM^w)k + D , . 

\ w= -Ki, sin + LM^w)k - ^M^w + {M^ s ^'^'^^ 

and taking away again from the first row and M'^ from the second row 
gives 

{ u= - {AA^y^A ic sin fc(v4^M + LM^w)k + D 



w= -{MM^y^MicSiukiA^u + LM^w)k - + j;S 



(34) 
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In this form it is clear that external currents D act as a bias only for the phase 
sector while external magnetic fluxes s act as a bias only for the mesh current 
sector. This is consistent with the physics of the system. Phases can grow 
indefinitely for suitable current biases and cannot be stopped by magnetic 
counterfields. Mesh currents instead are subject to a confining potential and 
an external bias cannot let them grow indefinitely, it can just modify their 
confined dynamics. In this case, it is seen from Eq. (|3^) that only an external 
magnetic field and not a bias current can change the position of mesh current 
minima, which are the flux dance steps. 

6 Conclusions 

The first conclusion is that in a mesh analysis a Josephson junction array 
SIMF model is described by an equation like Eq.(|T8|). For a circuit like b) 
the interaction range of phases is infinite so that it may be not necessary to 
take into account mutual inductance effects to explain finer system features 
like those studied in Ref. 0] or in Ref. [0]. The second conclusion is that the 
L —>■ limit of the SIMF model is actually the XY model and with Eq.(|30|) 
this limit can be studied analytically in a way easier than thought before. 
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